
%% Saturn
warning off
clc;clear;
format longg
options_negx = odeset('RelTol',1e-13,'AbsTol',1e-13,'Events',@NegXcrossing);
options=odeset('RelTol',1e-13,'AbsTol',1e-13);
options_fmincon = optimoptions('fmincon', 'MaxFunctionEvaluations', 5000, ...
    'MaxIterations', 5000,'ConstraintTolerance',1e-6,'Algorithm','Interior-point','stepTolerance',1e-21);
%% Inputs.
G           = 6.674e-11 * ((1/1000)^3); % Gravitational parameters

mass_central = 1.989e30;
%[] Mass of central planet

mass_moon = 5.972e24;         %Callisto
%[] mass of moons

DU = 151.73e6;           %Callisto
%[km] Semi-major axis of Moons, used as distance units for each CRTBP
%environment.

moonName = {'Luna'};
%[] Name of the Moons

thetao = 0;
%[] Initial position of the moon relative to the first point of aries

GM_central = mass_central*G;
%[] Gravitational parameters

GM_moon = mass_moon*G;
%[] Gravitational parameters

N = length(mass_moon);
%[] Number of Moons

u = zeros(N,1);
for ii = 1:N
    u(ii) = GM_moon(ii)/(GM_moon(ii)+GM_central);
end
%[] Gravataional ratio

TU = zeros(N,1);
VU = zeros(N,1);
for ii = 1:N
    TU(ii) = sqrt(DU(ii)^3/GM_central);
    VU(ii) = DU(ii)/TU(ii);
end
%[] Time constant for each crtbp system

theta_dot = zeros(1,N);
for ii = 1:length(theta_dot)
    theta_dot(ii) = sqrt(GM_central*DU(ii))/DU(ii)^2;
end
%[] Theta dot for each planet.


planetNumb = 1;
[L1,L2,L3,L4,L5] = librationPoints(u(planetNumb));
L_ = [L1,L2,L3,L4,L5];
for ii = 1:length(L_)
    L_pos = L_(:,ii);
    J_L(ii) = jacobiConst(L_pos,zeros(3,1),u(planetNumb));
end

%%
InclinationRotationAngle = 7.25;
RAANRotationAngle = 73.67+0.014*(2023-1850);
PoleDirection = Rotation([0;0;1],RAANRotationAngle,'Degrees')*Rotation([1;0;0],InclinationRotationAngle,'Degrees')*[0;0;1];
dt = 0.01;
startphase=2*pi*(4.5/12); % changing the starting position of the periodic 
              % orbit to show changes in the pole visibility based on the Satellite's position wrt the Sun's pole
              %%%% Optimal for L5 is 2*pi*(4.5/12)
              
%%
load('SE_L5_Vertical.mat');
inc_ind=1;
theta_list = [85,80,75,70,65,50];
% figure;hold on
for ii = 1:2:length(IC_L5_VL)
    IC = IC_L5_VL(:,ii);
    T = T_L5_VL(ii);
    [t,S] = ode113(@(t,S)CR3BP_n(t,S,u),[0:dt:T],IC,options);
    S = S';

    S_i = zeros(size(S));
    for ii = 1:length(S)
        S_i(:,ii) = C2I_primary(S(:,ii),u,DU,VU,t(ii)+startphase);
    end
    COE_temp = State2Coe(S_i(:,1),GM_central);
    inc_info(inc_ind) = COE_temp(3);
    theta_ind=1;
    for k = 1:length(theta_list)
        for ii = 1:length(S)
            angle = acosd(dot(PoleDirection,S_i(1:3,ii))/(norm(PoleDirection)*norm(S_i(:,ii))));
            if angle < theta_list(k)
                visible(ii) = 1;
            else
                visible(ii) = 0;
            end
        end
        VisibleTime(inc_ind,theta_ind) = sum(visible)*dt*TU/3600/24;
        if VisibleTime(inc_ind,theta_ind) == 0
            VisibleTime(inc_ind,theta_ind) = nan;
        end
        theta_ind=theta_ind+1;
    end
    inc_ind=inc_ind+1;
    clear S_i S
end




figure;
hold on

for ii = 1:length(theta_list)
plot(inc_info,VisibleTime(:,ii),'o-','MarkerIndices',[1:size(visible,2)]);
end
xlabel('Inclination')
ylabel('Pole access time [Days]');
set(gca,'FontSize',16)
legend('\theta=85','\theta=80','\theta=75','\theta=70','\theta=65')